FINAL PROJECT

#1. install necessary libraries

library(Seurat)
library(dplyr)
library(ggplot2)
library(patchwork)

options(future.globals.maxSize = 3 * 1024^3)  # Set to 3 GiB

#2 Load Control 10x Data

control = Load10X_Spatial(
     "/Users/cougerjaramillo/Desktop/MS/Intersession 2026/Applications/counts_and_images/1-1",
     filename = "filtered_feature_bc_matrix.h5",
     assay = "Spatial",
     slice = "slice1",
     bin.size = NULL,
     filter.matrix = TRUE,
     to.upper = FALSE,
     image = NULL
)
# Set orig.ident of control to 0
control@meta.data$orig.ident <- "0"
SpatialFeaturePlot(control, features = "nCount_Spatial") + theme(legend.position = "right")

control

An object of class Seurat 
36601 features across 3742 samples within 1 assay 
Active assay: Spatial (36601 features, 0 variable features)
 1 layer present: counts
 1 spatial field of view present: slice1

AD = Load10X_Spatial(
     "/Users/cougerjaramillo/Desktop/MS/Intersession 2026/Applications/counts_and_images/2-8",
     filename = "filtered_feature_bc_matrix.h5",
     assay = "Spatial",
     slice = "slice2",
     bin.size = NULL,
     filter.matrix = TRUE,
     to.upper = FALSE,
     image = NULL
 )
# Set orig.ident of AD to 1
AD@meta.data$orig.ident <- "1"
SpatialFeaturePlot(AD, features = "nCount_Spatial") + theme(legend.position = "right")

AD

An object of class Seurat 
36601 features across 3445 samples within 1 assay 
Active assay: Spatial (36601 features, 0 variable features)
 1 layer present: counts
 1 spatial field of view present: slice2

#3. Pre-processing

#control
#Mitochondrial DNA counts
control[["percent.mt"]] <- PercentageFeatureSet(control, pattern = "^MT-")#
#Visualize
VlnPlot(control, features = c("nFeature_Spatial", "nCount_Spatial", "percent.mt"), ncol = 3)

plot1 <- FeatureScatter(control, feature1 = "nCount_Spatial", feature2 = "percent.mt")
plot2 <- FeatureScatter(control, feature1 = "nCount_Spatial", feature2 = "nFeature_Spatial")
plot1

plot2

#Subset for high quality data
control <- subset(control, subset = nFeature_Spatial > 500 & nFeature_Spatial < 6000 & percent.mt < 40)
#Log normalize
control <- NormalizeData(control, normalization.method = "LogNormalize", scale.factor = 10000)

Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

#Feature Selection
control <- FindVariableFeatures(control, selection.method = "vst", nfeatures = 2000)

Calculating gene variances
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Calculating feature variances of standardized and clipped values
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

# Identify the 10 most highly variable genes
top10 <- head(VariableFeatures(control), 10)
# plot variable features with and without labels
plot1 <- VariableFeaturePlot(control)
plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE)
plot1

plot2

#Scale Data
all.genes <- rownames(control)
control <- ScaleData(control, features = all.genes)

  |                                          
  |                                    |   0%
  |                                          
  |=                                   |   3%
  |                                          
  |==                                  |   5%
  |                                          
  |===                                 |   8%
  |                                          
  |====                                |  11%
  |                                          
  |=====                               |  14%
  |                                          
  |======                              |  16%
  |                                          
  |=======                             |  19%
  |                                          
  |========                            |  22%
  |                                          
  |=========                           |  24%
  |                                          
  |==========                          |  27%
  |                                          
  |===========                         |  30%
  |                                          
  |============                        |  32%
  |                                          
  |=============                       |  35%
  |                                          
  |==============                      |  38%
  |                                          
  |===============                     |  41%
  |                                          
  |================                    |  43%
  |                                          
  |=================                   |  46%
  |                                          
  |==================                  |  49%
  |                                          
  |==================                  |  51%
  |                                          
  |===================                 |  54%
  |                                          
  |====================                |  57%
  |                                          
  |=====================               |  59%
  |                                          
  |======================              |  62%
  |                                          
  |=======================             |  65%
  |                                          
  |========================            |  68%
  |                                          
  |=========================           |  70%
  |                                          
  |==========================          |  73%
  |                                          
  |===========================         |  76%
  |                                          
  |============================        |  78%
  |                                          
  |=============================       |  81%
  |                                          
  |==============================      |  84%
  |                                          
  |===============================     |  86%
  |                                          
  |================================    |  89%
  |                                          
  |=================================   |  92%
  |                                          
  |==================================  |  95%
  |                                          
  |=================================== |  97%
  |                                          
  |====================================| 100%

#check final object for QC
control

An object of class Seurat 
36601 features across 3657 samples within 1 assay 
Active assay: Spatial (36601 features, 2000 variable features)
 3 layers present: counts, data, scale.data
 1 spatial field of view present: slice1

#AD
#Mitochondrial DNA counts
AD[["percent.mt"]] <- PercentageFeatureSet(AD, pattern = "^MT-")#
#Visualize
VlnPlot(AD, features = c("nFeature_Spatial", "nCount_Spatial", "percent.mt"), ncol = 3)

plot1 <- FeatureScatter(AD, feature1 = "nCount_Spatial", feature2 = "percent.mt")
plot2 <- FeatureScatter(AD, feature1 = "nCount_Spatial", feature2 = "nFeature_Spatial")
plot1

plot2

#Subset for high quality data
AD <- subset(AD, subset = nFeature_Spatial > 500 & nFeature_Spatial < 6000 & percent.mt < 40)
#Log normalize
AD <- NormalizeData(AD, normalization.method = "LogNormalize", scale.factor = 10000)

Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

#Feature Selection
AD <- FindVariableFeatures(AD, selection.method = "vst", nfeatures = 2000)

Calculating gene variances
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Calculating feature variances of standardized and clipped values
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

# Identify the 10 most highly variable genes
top10 <- head(VariableFeatures(AD), 10)
# plot variable features with and without labels
plot1 <- VariableFeaturePlot(AD)
plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE)
plot1

plot2

#Scale Data
all.genes <- rownames(AD)
AD <- ScaleData(AD, features = all.genes)

  |                                          
  |                                    |   0%
  |                                          
  |=                                   |   3%
  |                                          
  |==                                  |   5%
  |                                          
  |===                                 |   8%
  |                                          
  |====                                |  11%
  |                                          
  |=====                               |  14%
  |                                          
  |======                              |  16%
  |                                          
  |=======                             |  19%
  |                                          
  |========                            |  22%
  |                                          
  |=========                           |  24%
  |                                          
  |==========                          |  27%
  |                                          
  |===========                         |  30%
  |                                          
  |============                        |  32%
  |                                          
  |=============                       |  35%
  |                                          
  |==============                      |  38%
  |                                          
  |===============                     |  41%
  |                                          
  |================                    |  43%
  |                                          
  |=================                   |  46%
  |                                          
  |==================                  |  49%
  |                                          
  |==================                  |  51%
  |                                          
  |===================                 |  54%
  |                                          
  |====================                |  57%
  |                                          
  |=====================               |  59%
  |                                          
  |======================              |  62%
  |                                          
  |=======================             |  65%
  |                                          
  |========================            |  68%
  |                                          
  |=========================           |  70%
  |                                          
  |==========================          |  73%
  |                                          
  |===========================         |  76%
  |                                          
  |============================        |  78%
  |                                          
  |=============================       |  81%
  |                                          
  |==============================      |  84%
  |                                          
  |===============================     |  86%
  |                                          
  |================================    |  89%
  |                                          
  |=================================   |  92%
  |                                          
  |==================================  |  95%
  |                                          
  |=================================== |  97%
  |                                          
  |====================================| 100%

#check final object for QC
AD

An object of class Seurat 
36601 features across 3366 samples within 1 assay 
Active assay: Spatial (36601 features, 2000 variable features)
 3 layers present: counts, data, scale.data
 1 spatial field of view present: slice2

#4. Merge the control & AD

#Merging Slices
brain.merge <- merge(control, AD, project = "brain")
#VariableFeatures(brain.merge) <- c(VariableFeatures(control), VariableFeatures(AD))
#Re-Scale Data
all.genes <- rownames(brain.merge)
brain.merge <- ScaleData(brain.merge, features = all.genes)

  |                                          
  |                                    |   0%
  |                                          
  |=                                   |   3%
  |                                          
  |==                                  |   5%
  |                                          
  |===                                 |   8%
  |                                          
  |====                                |  11%
  |                                          
  |=====                               |  14%
  |                                          
  |======                              |  16%
  |                                          
  |=======                             |  19%
  |                                          
  |========                            |  22%
  |                                          
  |=========                           |  24%
  |                                          
  |==========                          |  27%
  |                                          
  |===========                         |  30%
  |                                          
  |============                        |  32%
  |                                          
  |=============                       |  35%
  |                                          
  |==============                      |  38%
  |                                          
  |===============                     |  41%
  |                                          
  |================                    |  43%
  |                                          
  |=================                   |  46%
  |                                          
  |==================                  |  49%
  |                                          
  |==================                  |  51%
  |                                          
  |===================                 |  54%
  |                                          
  |====================                |  57%
  |                                          
  |=====================               |  59%
  |                                          
  |======================              |  62%
  |                                          
  |=======================             |  65%
  |                                          
  |========================            |  68%
  |                                          
  |=========================           |  70%
  |                                          
  |==========================          |  73%
  |                                          
  |===========================         |  76%
  |                                          
  |============================        |  78%
  |                                          
  |=============================       |  81%
  |                                          
  |==============================      |  84%
  |                                          
  |===============================     |  86%
  |                                          
  |================================    |  89%
  |                                          
  |=================================   |  92%
  |                                          
  |==================================  |  95%
  |                                          
  |=================================== |  97%
  |                                          
  |====================================| 100%

#Visualize Results
plot1 <- VlnPlot(brain.merge, features = "nCount_Spatial", pt.size = 0.1) + NoLegend()
plot2 <- SpatialFeaturePlot(brain.merge, features = "nCount_Spatial") + theme(legend.position = "right")
plot1

plot2

#5 # Dimensionality reduction, clustering, and visualization

#PCA
brain.merge <- RunPCA(brain.merge, assay = "Spatial", verbose = FALSE)

Error: vector memory limit of 16.0 Gb reached, see mem.maxVSize()

#6 Identifying Cortical Layers by Heatmpap of canonical markers

# Define the list of known markers organized by layer
known_markers_list <- list(
    Layer1 = c("SPARC", "MALAT1", "CXCL14", "ADIRF"),
    Layer2 = c("ENC1", "HPCAL1", "CALB2"),
    Layer3 = c("HOPX"),
    Layer4 = c("RORB", "NUAK1"),
    Layer5 = c("TUBA1B", "TMSB10", "SYT1", "STMN2", "MAP1B", "SNAP25", "PCP4", "UCHL1", "CLSTN2", "SNCG", "SMYD2", "RTN1", "STMN1", "RTN3", "NEFL", "SNRPN", "BASP1", "SYN2", "CALM3", "ENO2", "SNCA", "GAP43", "NAPB", "ELAVL4", "FXYD6", "NDRG4", "NRSN1", "RAB3C", "UHMK1", "SNAP91", "ATP6V1B2", "SLC25A22"),
    Layer6 = c("MAP2K1", "DIRAS2", "KRT17"),
    LayerWM = c("MBP", "MOBP", "CLDND1", "BCAS1", "MTURN", "PAQR6", "HIPK2", "DYNC1LI2") 
)
# Combine all markers into a single vector and create a corresponding layer vector
combined_markers <- unlist(known_markers_list)
layer_labels <- rep(names(known_markers_list), times = sapply(known_markers_list, length))
# Ensure combined markers exist in the Seurat object
valid_markers <- combined_markers[combined_markers %in% rownames(brain.merge)]
valid_layer_labels <- layer_labels[combined_markers %in% rownames(brain.merge)]
# Check if we have valid markers before plotting
if (length(valid_markers) > 0) {
    # Create a named vector for y-axis annotations
    names(valid_layer_labels) <- valid_markers
    
    # Generate the heatmap
    heatmap_plot <- DoHeatmap(brain.merge, features = valid_markers) + 
        NoLegend() + 
        scale_y_discrete(labels = valid_layer_labels) +  # Update y-axis labels
        ggtitle("Combined Heatmap of Known Markers by Layer")  # Title for the heatmap
    
    # Display the heatmap
    print(heatmap_plot)
} else {
    cat("No valid markers found for heatmap generation.\n")
}

#7. Subset the Data by Region

#Subset out anatomical regions
## Subset to clusters of interest
WM <- subset(brain.merge, idents = c(3, 7))

# Extract and attach spatial coordinates from image 1
centroids <- WM[["slice1"]]@boundaries$centroids
coords <- setNames(as.data.frame(centroids@coords), c("x", "y"))
rownames(coords) <- centroids@cells
WM$x <- coords[colnames(WM), "x"]
WM$y <- coords[colnames(WM), "y"]

# After subsetting, we renormalize WM
WN <- NormalizeData(WM, normalization.method = "LogNormalize", scale.factor = 10000)

Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

WM <- RunPCA(WM, verbose = FALSE)
SpatialDimPlot(WM, label = TRUE, label.size = 3)

#8 WM DEGs

# Join Layers
WM <- JoinLayers(WM)

# Set the identity classes based on the 'orig.ident' column
WM <- SetIdent(WM, value = WM$orig.ident)


# Compare cells with orig.ident values of 0 and 1
WM.AD.markers <- FindMarkers(WM, ident.1 = 0, ident.2 = 1)

  |                                                  | 0 % ~calculating  
  |+                                                 | 1 % ~01m 07s      
  |+                                                 | 2 % ~01m 06s      
  |++                                                | 3 % ~01m 05s      
  |++                                                | 4 % ~01m 04s      
  |+++                                               | 5 % ~01m 03s      
  |+++                                               | 6 % ~01m 02s      
  |++++                                              | 7 % ~01m 01s      
  |++++                                              | 8 % ~59s          
  |+++++                                             | 9 % ~58s          
  |+++++                                             | 10% ~57s          
  |++++++                                            | 11% ~56s          
  |++++++                                            | 12% ~55s          
  |+++++++                                           | 13% ~54s          
  |+++++++                                           | 14% ~54s          
  |++++++++                                          | 15% ~53s          
  |++++++++                                          | 16% ~52s          
  |+++++++++                                         | 17% ~52s          
  |+++++++++                                         | 18% ~51s          
  |++++++++++                                        | 19% ~50s          
  |++++++++++                                        | 20% ~50s          
  |+++++++++++                                       | 21% ~49s          
  |+++++++++++                                       | 22% ~49s          
  |++++++++++++                                      | 23% ~48s          
  |++++++++++++                                      | 24% ~47s          
  |+++++++++++++                                     | 25% ~47s          
  |+++++++++++++                                     | 26% ~46s          
  |++++++++++++++                                    | 27% ~46s          
  |++++++++++++++                                    | 28% ~45s          
  |+++++++++++++++                                   | 29% ~44s          
  |+++++++++++++++                                   | 30% ~44s          
  |++++++++++++++++                                  | 31% ~44s          
  |++++++++++++++++                                  | 32% ~43s          
  |+++++++++++++++++                                 | 33% ~42s          
  |+++++++++++++++++                                 | 34% ~42s          
  |++++++++++++++++++                                | 35% ~41s          
  |++++++++++++++++++                                | 36% ~41s          
  |+++++++++++++++++++                               | 37% ~40s          
  |+++++++++++++++++++                               | 38% ~39s          
  |++++++++++++++++++++                              | 39% ~39s          
  |++++++++++++++++++++                              | 40% ~38s          
  |+++++++++++++++++++++                             | 41% ~37s          
  |+++++++++++++++++++++                             | 42% ~37s          
  |++++++++++++++++++++++                            | 43% ~36s          
  |++++++++++++++++++++++                            | 44% ~35s          
  |+++++++++++++++++++++++                           | 45% ~35s          
  |+++++++++++++++++++++++                           | 46% ~34s          
  |++++++++++++++++++++++++                          | 47% ~33s          
  |++++++++++++++++++++++++                          | 48% ~33s          
  |+++++++++++++++++++++++++                         | 49% ~32s          
  |+++++++++++++++++++++++++                         | 50% ~31s          
  |++++++++++++++++++++++++++                        | 51% ~31s          
  |++++++++++++++++++++++++++                        | 52% ~30s          
  |+++++++++++++++++++++++++++                       | 53% ~30s          
  |+++++++++++++++++++++++++++                       | 54% ~29s          
  |++++++++++++++++++++++++++++                      | 55% ~28s          
  |++++++++++++++++++++++++++++                      | 56% ~28s          
  |+++++++++++++++++++++++++++++                     | 57% ~27s          
  |+++++++++++++++++++++++++++++                     | 58% ~27s          
  |++++++++++++++++++++++++++++++                    | 59% ~26s          
  |++++++++++++++++++++++++++++++                    | 60% ~25s          
  |+++++++++++++++++++++++++++++++                   | 61% ~25s          
  |+++++++++++++++++++++++++++++++                   | 62% ~24s          
  |++++++++++++++++++++++++++++++++                  | 63% ~23s          
  |++++++++++++++++++++++++++++++++                  | 64% ~23s          
  |+++++++++++++++++++++++++++++++++                 | 65% ~22s          
  |+++++++++++++++++++++++++++++++++                 | 66% ~21s          
  |++++++++++++++++++++++++++++++++++                | 67% ~21s          
  |++++++++++++++++++++++++++++++++++                | 68% ~20s          
  |+++++++++++++++++++++++++++++++++++               | 69% ~19s          
  |+++++++++++++++++++++++++++++++++++               | 70% ~19s          
  |++++++++++++++++++++++++++++++++++++              | 71% ~18s          
  |++++++++++++++++++++++++++++++++++++              | 72% ~17s          
  |+++++++++++++++++++++++++++++++++++++             | 73% ~17s          
  |+++++++++++++++++++++++++++++++++++++             | 74% ~16s          
  |++++++++++++++++++++++++++++++++++++++            | 75% ~15s          
  |++++++++++++++++++++++++++++++++++++++            | 76% ~15s          
  |+++++++++++++++++++++++++++++++++++++++           | 77% ~14s          
  |+++++++++++++++++++++++++++++++++++++++           | 78% ~14s          
  |++++++++++++++++++++++++++++++++++++++++          | 79% ~13s          
  |++++++++++++++++++++++++++++++++++++++++          | 80% ~12s          
  |+++++++++++++++++++++++++++++++++++++++++         | 81% ~12s          
  |+++++++++++++++++++++++++++++++++++++++++         | 82% ~11s          
  |++++++++++++++++++++++++++++++++++++++++++        | 83% ~11s          
  |++++++++++++++++++++++++++++++++++++++++++        | 84% ~10s          
  |+++++++++++++++++++++++++++++++++++++++++++       | 85% ~09s          
  |+++++++++++++++++++++++++++++++++++++++++++       | 86% ~09s          
  |++++++++++++++++++++++++++++++++++++++++++++      | 87% ~08s          
  |++++++++++++++++++++++++++++++++++++++++++++      | 88% ~07s          
  |+++++++++++++++++++++++++++++++++++++++++++++     | 89% ~07s          
  |+++++++++++++++++++++++++++++++++++++++++++++     | 90% ~06s          
  |++++++++++++++++++++++++++++++++++++++++++++++    | 91% ~06s          
  |++++++++++++++++++++++++++++++++++++++++++++++    | 92% ~05s          
  |+++++++++++++++++++++++++++++++++++++++++++++++   | 93% ~04s          
  |+++++++++++++++++++++++++++++++++++++++++++++++   | 94% ~04s          
  |++++++++++++++++++++++++++++++++++++++++++++++++  | 95% ~03s          
  |++++++++++++++++++++++++++++++++++++++++++++++++  | 96% ~02s          
  |+++++++++++++++++++++++++++++++++++++++++++++++++ | 97% ~02s          
  |+++++++++++++++++++++++++++++++++++++++++++++++++ | 98% ~01s          
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 99% ~01s          
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=01m 02s

# View the results
head(WM.AD.markers)